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Abstract 

In close analogy to diffusion limited aggregation (DLA) and in- 
spired by a work of Roux, a random walker algorithm is constructed 
to solve the problem of crack growth in an elastic medium. In contrast 
to conventional lattice approaches, the stress field is not calculated 
throughout the whole medium, but random walkers are used to detect 
only the hot sites on the surface of the crack. There, an analytically 
calculated Green's function is used to determine the stress field. The 
complicated boundary condition on the crack surface is simulated by a 
special sticking-rule walk. Using this new method we generate crack- 
clusters up to sizes of 20,000 particles on simple workstations within 
reasonable time. We simulate several different boundary conditions, 
like uniaxial tensile, pure shear and isotropic tensile load. We study 
the influence of several parameters representing the material strength 
or fatigue. Furthermore we study the effect of anisotropic walks. As a 
result, we reproduce with this new model typical experimental crack 
shapes and are able to simulate the essential features of realistic cracks. 



"HLRZ preprint 106/92 



1 



1 Introduction 



The study of non-equilibrium growth models has been extremely popular in 
the past few years. Especially the Laplacian growth phenomenon has been 
studied extensively and applications of it can be found in a large variety 
of fields, e. g. electrodeposition |], fluid-fluid displacement @ or growth 
of bacteria colonies ||. One reason for its popularity is, that a numerical 
simulation of this process can be done extremely efficiently using diffusion- 
limited aggregation (DLA) ||. Using highly optimized programs || |J it 
is possible to grow huge DLA clusters containing up to 50 million particles 

10. B- 

However, this paper is not devoted to DLA, but inspired by this efficiency 
and success of DLA and by a work of Roux || we will attack the problem of 



crack growth in an elastic medium |10], |TTJ, [12]]. Here, to an elastic material 
a certain load is applied and eventually cracks develop, which grow until 
the material breaks apart. Although highly developed numerical techniques, 
which are mainly based on lattice models, exist 0, [14]], the simulation of 
crack growth is even on today's supercomputers restricted to a few thousand 
broken bonds or even less. Thus, new techniques are urgently needed, which 
might allow for large crack growth simulations maybe on workstations. 

As was pointed out by many researchers |15|, [L6|] , there exist some remark- 
able formal similarities between the vectorial elasticity problem — which is 
governed by the Lame equation — and the scalar electricity problem — which 
is governed by the Laplace equation. It was even shown by Roux, that it is 
formally possible to construct a Green's function of the Lame equation using 
a certain type of vectorial random walkers in a similar way one constructs 
a Green's function of the Laplace equation. The latter one can of course 
simply be interpreted as a density of random walkers and is thus much easier 
to implement than a formal, vectorial random walker, which is much more 
difficult to interpret. However, although Roux's ideas are very clear, an ac- 
tual numerical implementation of them and a check of their practicability is 
still missing. 

In the current paper we will discuss his ideas and will present a new 
algorithm for crack growth. For this purpose we will construct a DLA-like 
random walker method, in which the random walkers are used to find the 
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"hot" (highly stressed) sites along the surface of a crack, i. e. those which 
are most likely to break. Here, 'DLA-like' is not meant in the sense of 
'probabilistic' used in [T7|, [TB], [2tJ to characterize these crack growth 



models, but in a sense that we explicitly use random walkers to transport 
information from the boundaries towards a crack. 

The use of random walkers will introduce disorder and noise into the 
simulation. But since real material are typically spatially disordered and 
inhomogeneous, this noise can be justified and interpreted physically as an 
inherit randomness. It does not have to be introduced by hand like in typical 
lattice models, but is an implicit feature of the model. 

The outline of this paper is as follows. In section ^| we will give a short 
introduction into the problem of crack growth. We will briefly present and 
extend Roux's ideas and we will discuss the reasons, why a direct numerical 
simulation of his method is not practical. In section [3] we will introduce a 
new model for crack growth by drawing an analogy to DLA. We will use 
random walkers and an analytically calculated Green's function to transport 
the information about surface forces towards a crack, which is represented 
as a cluster of particles. We will also discuss the treatment of the boundary 
condition along the crack, which turns out to be the most difficult ingredient 
of the model. Then, in section [|, we shall present the results of this model 
and will study several variations of it in order to understand the influence of 
the numerical parameters. Here, we will compare our model to experimental 
and analytical results. In section |5| we will critically discuss the method and, 
finally, in section ^ we will summarize our results. 



2 Relation between crack growth and ran- 
dom walks 

Consider a d dimensional infinite elastic medium with Poisson ratio v . In 
equilibrium this system can be described in terms of the stress field a which 

obeys 

Div^ = 0. (1) 
Together with Hooke's law this can be reformulated in terms of the dis- 
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placement field u and one obtains the Lame equation 



(1 - u{d - 1))V 2 u + V(V- u) = 0. 



(2) 



A crack in such a system can be interpreted as an additional force free 
surface, which thus has to obey the boundary condition 



Here one has to point out, that unlike in DLA this boundary condition 
at the crack surface is not formulated in terms of the displacement field - 
which corresponds to the potential in the Laplace equation — but in terms 
of the stress field, which is related via Hooke's law to a spatial derivative of 
u . As we will show later, this peculiarity requires a special treatment of the 
crack surfaces. 

Apart from those equilibrium conditions one can formulate the condition 
for the growth of a crack in terms of the stress component parallel to the 
crack surface a\\ . This has to exceed a critical value a c , which is directly 
connected to the intermolecular cohesion forces, and then the growth velocity 
of the crack v n is determined by 



where r] is a heuristic parameter which is often simply set to one. This 
rule for crack growth only takes into account crack growth through cleavage 
and completely ignores bending terms. The equations (||, [|, |) formulate the 
problem of crack growth as a nonlinear moving boundary problem. 

To give a closed picture of the problem we shall now recall very briefly 
the main ideas of Roux || in order to be able to extend and comment them. 
The fundamental idea of his work is to define a Markovian random walk 
process that constructs a matrix, which finally leads to a Green's function of 
the Lame equation. This is done in complete analogy to the Laplacian case, 
which we will find as a special case of the following derivation. 

Therefore, consider an elastic material with an unknown displacement 
field u , which is a result of the boundary conditions. Our task is to obtain 
this field at point x in terms of the boundary displacements. Therefore 



a- n = 0. 



(3) 




(4) 
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consider the following Markovian process: A random walker starts at point 
x and makes a first, elementary jump of length r into a random direction 
e 1 . After averaging over all possible directions e\ one defines at xa new 

vector (Q r u)( x) as 



(Q r u)(x) — (Q( ei)- u( x + r e^j ^ (5) 

where the elementary step matrix Q{ e) is defined as 

Q( e) . = (ck-1 + (1 - a)-d- e <g> e) = (a-^- + (1 - aj-d-aej) (6) 

and a is a parameter that will be determined later. 

As usual in Markov processes this elementary step is iterated. Conse- 
quently, after performing N independent steps - each of length r - one defines 
at x the new vector 



(Q? u)( x) = (Q( ejv )---£( ei )- u( x + r Cl + ■ ■ ■ + r e N )) ^ . (7) 

By considering its spatial Fourier transform 

(Q* u ){k)= J (Q? u) ( x) exp (-i k x)d d x = k)- u{ k) (8) 

one can evaluate the limit of vanishing step length r— >0 and diverging 
number of steps N^oo , while keeping t = Nr 2 fixed. One easily obtains 
with the central-limit theorem 

lim ( k) = Q t (k) = exp (- 2d ^ +2 ^ [i 2a + d ) k \ + 2d ( l -a)k® k 



r—>0 



= exp (m) • 

(9) 

After transforming this quantity back into real space, it is easily seen that 
the Fourier transform 

(Q* u)(x) = J Q*(r)- u(x+ r)d d r (10) 
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provides a solution of the following general equation for the vector field 
V(r,t) 



d V _ (2a + d)V 2 V + 2d(l - q)V(V- V) 
dt ~ 2d(d + 2) { } 

The stationary limit d V/dt — of this equation obviously provides the 
Lame equation if one chooses a appropriately 

d(l-2u(d-l)) 
a 2{d + 1 - dv{d - 1)) { ' 

and identifies V( r,t) with the displacement field u( r,t) . 

In summary, in this formulation the problem of determining a time de- 
pendent Green's function Q t ( r) , which is used in (|TUD, has been formulated 
in terms of a random walk. Obviously eq. (|7|) shows the type of Markov 
process that one has to perform in order to solve Lame's equation. Anyhow, 
we shall show in the following that although exact, it is not practical to 
construct an algorithm directly from (0). 

A direct and naive implementation of ([7p would involve the following 
steps. One would launch a random walker at a boundary site r j = x + r at 
which the displacement u is known. The walker would undergo an isotropic 
random walk until it touches a crack at point r/ = x . During the walk one 
would calculate the product of the elementary step matrices Q( e^r) • • • Q{ ei) 
and would finally average at x over the vectors Q( e^) ■ • • Q( ex)- u( r^) of all 
incoming random walkers. One would then have to extrapolate to the limit 
r— >0 — which at least in this algorithm necessarily leads to N^oo — and 
according to (|10D this procedure should finally yield the required solution. 

But although this process strongly reminds one of the process to construct 
the Green's function of the Laplace equation with random walkers, there 
are several differences which result in important difficulties. Unlike in DLA 
the history of each walk is here extremely important. The product matrix 
is explicitly dependent on each individual step e^ and therefore the path 
on which the walker travels from its initial to its final point is extremely 
important. Thus, the contribution of each individual walk to the average (0) 
is arbitrarily small. 
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On the other hand it is necessary to take the limit iV— >oo and thus 
the number of factors to calculate the product matrix grows rapidly. But 
since the eigenvalues of Q( e) are not equal to unity the product matrix 
Q{ ejv) • • • Q( ei) has an uncontrollable norm, which either diverges or van- 
ishes exponentially fast. From this follows, that also the modulus of the re- 
sulting vector Q{ e^) • • -Q( ei)- u( rj) is uncontrollable. Similar arguments 
show immediately that not only the modulus, but also the direction of the 
resulting vector is uncontrollable. The consequence of these considerations 
is, that a naive implementation of the algorithm suggested by (0) results in 
enormous statistical fluctuations of both the modulus and the direction of 
(^Qr uj ( x) , which cannot be suppressed in numerical simulations. 

Another important point is, that the previously described algorithm ne- 
glects the presence of a crack, which substantially changes the stress field. It 
is also not easily possible to include the boundary condition at the crack (|3|), 
because it is expressed in terms of the stress field, which is not included in the 
previous derivation. Of course one could reformulate the boundary condition 
in terms of the displacement field, but this requires spatial derivatives of the 
displacement field which are even less accessible than the displacement itself. 

However, the previous considerations showed that one purpose of the 
random walkers is to construct the matrix Q t ( r) , which is the basis of fliPf) . 
But in an infinite system this matrix, which is actually already the average 
of all product matrices of walks starting at the origin and terminating at r 
, can be calculated analytically. Thus, one could evaluate ( JT0|) directly using 
random walks, which would then only be used to determine the initial and 
final positions x and x + r . Here, the history of each individual walk is no 
longer relevant and thus one does not obtain the huge fluctuations described 
before. 

Q t ( r) is obtained by a Fourier transformation of (§) into real space (see 
Appendix). Finally one finds 

2*( r) = hf s (r,t) + f a (r,t)}-l - J-/ a ( r ,f)-2-^ (13) 
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where 



exp (-r 2 /4Ait) exp (-r 2 /4A 2 t) 
Mr,tj " AX 1 t + AX 2 t 



fm{r , t)= fi + ^) °P(-fW - ( 1 + ^ Vx P ( ,74A 2 t) ^ 
\ H / 4Ait \ r 2 / 4A 2 t 

This result has to be commented in several respects. 

It is clearly seen from the form of the two functions f s and /„ , that 
the underlying elementary process is a random walk. They consist of two 
Gaussians, which are the result of the diffusion of free particles launched at 
the origin of an infinite domain. The limiting case of the pure Laplacian is 
obtained by setting v— > — oo . In this case Q ( r) simplifies to the Green's 
function of the corresponding Laplace problem, which is a pure Gaussian. 

Since one requires t = Nr 2 to be fixed, one has introduced a new, free 
time variable. But none of the recent and efficient implementations of DLA 
algorithms uses an explicit time definition. In DLA the time it takes for a 
walker to reach a cluster is not relevant. Here, k explicitly enters Q l ( r) . 

The previously calculated Q l ( r) is only valid in an infinite domain. In 
(0) one averages essentially over all possible random walks that connect the 
origin with the point r . But in the presence of a crack, many of the walks 
terminate at the crack before reaching r . Also those walks are taken into 
account in eqs. (|7[-|T0l) to calculate Q\ r) . This deficiency could be altered 
by considering only those walks, which do not leave a restricted area, like 
in calculations of "first-passage-times", but this case cannot be calculated 
analytically. 

Equation (pi]) leads to the solution of the Lame equation only in the 
stationary limit t^oo . But in this limit Q l { r) vanishes identically, which 
means that one only obtains the trivial solution u{ r) = . 

As a result one obtains, that also a numerical simulation of (|l(]) — al- 
though technically possible — cannot be done, since mainly the time variable 
t cannot be interpreted in a simulation. 

The essential result of the previous considerations is, that the matrix 
Q\ r) plays the role of a Green's function of the time dependent problem (|Tll) 
and the original algorithm simply provides a prescription how to construct 
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it using random walks. On the other hand, one would in all cases end up 
with a Green's function in an infinite domain, which in fact can be calculated 
directly and analytically. 



3 Algorithm 

From the previous sections it is clear that the use of the original random 
walker method always ends up in the calculation of a Green's function for 
the displacement field. Both, the boundary condition at the crack surface 
and the breaking criterion, are usually given in terms of the stress field which 
is related to the spatial derivative of the displacement field using Hooke's law. 
Thus it seems natural to calculate the Green's function for the stress field 
analytically and use this to construct a random walker algorithm. 

Therefore, consider an infinite domain in the absence of cracks which 
is kept fixed at infinity. At the origin one applies a point force F which 
generates a stress field a*( r) throughout the whole system and is governed 



by Div a* = F-8( r) . Using the method of Kolossov-Mushkelishvili [21, 22 
we can calculate the stress field in the whole system and obtain 



*G) + *(? + *i) 



<y = m Q~^ + 4 ] (ir,) 



zc c 

-^ + x- 

z z z 



where we represent all 2d vectors as complex numbers 



x + zy 

(16) 



" 2tt(1 + X ) 
X = 3 - 4i/ 

and dt(z) ( ^s(z) ) denote the real (imaginary) part of z . Now one has to 
take into account that the probability for a random walker, that is launched 
from the center of a circle of radius r , to reach a certain element dip at the 
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perimeter of this circle is pdip = dip/2irr . This term has to be divided out 
of (15) and therefore we define 



g=2ir\z\%*. (17) 

Using this Green's function a, which is the basis of the following algo- 
rithm, we define a new method to determine the stress field in the medium: 

1. Since (|IJ) is linear, we can obtain the solution for a line of forces as 
a superposition of elementary point forces. Like in DLA, where each 
point which is kept at a nonzero potential is interpreted as a source of 
walkers, each point force acts also here as a source of random walkers. 
Random walkers are launched — like in DLA one at a time — with 
equal probabilities from each point where a force Fi is applied. All 
forces are usually of unit strength. Each walker undergoes an ordinary 
random walk until it touches a particle of the cluster located at point rf 
. Afterwards one calculates the stress tensor a( ry — Fi) according 
to flTTD and accumulates this stress in the counter of the appropriate 
surface element (see below). 

2. A crack is represented as an aggregate of particles of unit diameter. 
Like in DLA the growth of a crack has to be initiated by locating a 
seed particle at some point, which then acts as a single microcrack, from 
which the crack grows. Since we expect strong fluctuations during the 
growth of the crack we use a special scheme to reduce the fluctuations 
and get a better estimate of the stress field. It is well known from DLA 



that the introduction of anisotropic noise-reduction schemes [23 



immediately leads to strong lattice effects. Although similar schemes 



have not shown such effects in lattice models of crack growth [25[ we try 
to avoid possible problems, and therefore our scheme must not define 
any surface nor prefer any direction. The crack is represented by the 
particles as indicated in figure 1: Each particle represents two stress 
counters, one for each crack surface. Only the particles at the tip of 
a branch, i. e. those with only one parent particle but no children, 
have one single counter. In these counters the stresses a carried by the 
incoming walkers are accumulated. 
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3. Obviously the correct boundary condition (|) is not obeyed by simply 
calculating stresses according to fllTD , but an unbalanced force 

^unbal = £ • H ± 
remains. Therefore, we apply a virtual balancing force 

^bal = » ( 18 ) 

to the appropriate surface element at ry , which exactly compensates 
the spurious shear component r and stress component perpendicular 
to the crack surface a± . Thus it is in principle sufficient to accumulate 
at the crack only the stress component parallel to the crack surface an 
. But on the other hand we will show later that the shear component 
will be essential to determine the direction into which the crack grows 
and thus we accumulate at r/ the shear component as well as a\\ . As 
was said above, a random walker is started at each point where a force 
is applied. Consequently, the walker has to be restarted from r/ and 
is reinserted into the system, this time with the new force -P^al • This 
process of constantly reinserting the walker into the system is repeated 
until the modulus of i*bal drops below a certain threshold, which in 
our simulations is typically chosen as 10~ 8 . Afterwards the walker is 
discarded and a new walker is started. This repeated process can be 
interpreted as a method to relax the stress at the surface and can best 
be understood in terms of a boundary integral method [ 2!J . Physically 



this process is like a multiple scattering method. A similar situation 
occurs in DLA in simulations of viscous fingering with nonvanishing 



surface tension [27], [28], or in simulations of electrochemical deposi- 
tion with nonvanishing surface impedance . Both cases are mapped 
to a simulation of sticking probability DLA, in which the walker stops 
at the cluster with a probability smaller than unity. 

4. Next we have to define a growth rule. Here it has to be stressed that 
the growth of the crack must not be governed by the fact that the crack 
is touched by a walker, but only the accumulated stress tensor should 
determine the growth of the crack. Since the stress tensor is symmetric, 
it defines eigenvalues and eigenvectors, which determine the principal 
stresses and principal stress directions. The eigenvalues also allow to 
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distinguish between tensile and compressive stress, because both have 
a different sign. Now, we provide two material parameters: one for 
tension <j t and one for compression o c . They represent the material 
strength under the corresponding load. Especially the strength under 
tensile load a t can be related to the intermolecular cohesion forces. If 
either of them is exceeded by one eigenvalue, the crack grows and a 
new particle is added to the crack. The direction into which to put the 
new particle is determined by the eigenvector of the other eigenvalue. 
This breaking rule is chosen such that for a pure uniaxial tension the 
crack grows perpendicular to the direction of the force. In this growth 
rule only the principal stresses, which are purely tensile or compressive, 
are checked against the material strength. This situation is similar to 



the growth rule in the central-force model ||31|| , in which also no shear 
or bending modes are used to break a bond. 

5. Once the crack has grown it has to release its stress. Here, we study two 
different situations. In the first one the stress is only released locally, 
i. e. only the stress counter of the particle to which a new particle 
is added is cleared, while all other counters keep their values. This 
situation can physically be compared to a material in which fatigue 
plays an important role since even a small load can accumulate and 
break a sample if it is acting long enough. In the second situation 
we study a global relaxation. In this scheme we not only clear the 
counter of the particle at which the crack has grown, we also clear all 
counters in the entire cluster after the crack has grown by M particles. 
The parameter M is a fatigue parameter. The first situation can be 
expressed in terms of this parameter as M = oo . 



4 Results 

In the following we are going to present the results obtained by simulating 
the previous method. We simulated several different boundary conditions: 
uniaxial tension and compression, isotropic tension and pure shear (fig. ). In 
order to compare with experimental results we also simulated a four point 
shear geometry and another geometry used in experiments . We furthermore 
studied the influence of material strength and fatigue. 
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The length scale in all following simulations is the particle diameter and 
all applied forces have unit length, which defines the unit of the stress. 

4.1 Basic results 

In the following we will present the results of three different simulations with 
the boundary conditions uniaxial tension, isotropic tension and pure shear. 
The material strengths are typically chosen to be a t = 20 and o c = —20 (sc. 
the applied forces are of unit size). Thus, the system breaks equally well 
under compression and under tension. Since the stress components parallel 
to the crack surface of one walker are typically around 0.5 , at least N = AO 
walker have to touch a certain site before the cluster can grow. As the seed 
of all clusters we place two particles with unit diameter at the points (0,0.5) 
and (0,-0.5). This initial geometry is chosen arbitrarily, but simulations with 
other geometries show that the results depend neither on the number nor on 
the specific geometry of these seed particles. The Poisson ratio is v — 0.2 and 
the fatigue parameter is M = 10. All clusters have a size of 6, 000 particles. 

In figure 3 we show a typical result of a simulation of uniaxial tension. 
Here, forces are located along the lines f — (±2000, y), y G [—2000,2000] 
and point into the positive (negative) x direction. One obtains essentially a 
one-dimensional, straight crack without side branches that grows perpendic- 
ular to the direction of the force. This result is physical and can easily be 
observed in experimental and other numerical studies. 

In figure 4 we show a typical result of a simulation of pure shear. In this 
simulation all forces are again located along the lines r = (±2000, y). But 
now the forces are pointing into the positive (negative) y direction to produce 
a pure shear field. Although we use the same microcrack as in the previous 
simulation the shape of the crack is now drastically altered. The resulting 
cracks show a pronounced cross like shape. This main shape is a result of 
the pure shear field, whose stress tensor contains only off-diagonal elements, 
and thus the eigenvectors always point into a 45° angle from a surface. We 
have to remind, that this shape is not a lattice effect but a consequence of 
the applied forces. 

The reason for the symmetry of the cross-like shape is that the mate- 
rial strengths for tension and compression are symmetric, a c = —a t . The 
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case of asymmetric strength will be shown later. The cross-like shape under 
symmetric breaking conditions has been observed in other numerical studies 
[IIR III] on much smaller length scales. In experiments this shape is usually 



not observed, because the material strength for tension is usually smaller 
than for compression. 

One also observes small side branches which form right angles with the 
main branch. We have to stress, that these right angles cannot be the result 
of lattice effects because the whole method is formulated off-lattice. Our 
simulations are the first off-lattice simulations which show this behavior. The 
formation of right angles is completely determined by the cracking process 
and has a physical origin: Since there are no restoring forces acting on the 
crack surface, the boundary condition (|3|) only allows stresses parallel to the 
surface. Consequently side branches can only grow perpendicular to the main 
branches in the vicinity of the main crack. 

A typical result for a crack pattern under isotropic tension is shown in 
figure 5. Here, the radially outward pointing forces are located along a circle 
of radius 2000. One obtains a ramified structure which can clearly be dis- 
tinguished from DLA: The typical tip-splitting instability vanished and the 
side branches show again a tendency to grow perpendicular to their main 
branches. 



4.2 Fatigue 

As a first variation we want to study the influence of the fatigue parameter 
by setting M = oo while keeping all other parameters fixed. Physically 
we simulate now a material in which a small load can accumulate and can 
damage the material. Eventually this can lead to the formation of a large 
crack. 

Since we use M = oo we need much less walkers to reach the material 
strength and can simulate much larger cracks. Here, we show cracks contain- 
ing 20, 000 particles. 

Under uniaxial tension (fig. 6) we obtain again a narrow straight crack 
perpendicular to the applied load, but we observe the formation of many 
small side branches, which grow perpendicular to the main crack and which 
can themselves form perpendicular side branches. We have to stress again 
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that the formation of perpendicular side branches is not a result of lattice 
effects but quite physical. However, the formation of multiple side branches 
is not physical and thus we obtain, that the fatigue parameter is essential to 
describe the physical processes. 

Under shear load (fig. 7) we observe the formation of the cross-like main 
branches. Like in the previous case the introduction of fatigue leads also in 
this case to the formation of many small side branches which again may form 
side branches themselves. This leads eventually to a dendritic shape of the 
crack. 

Very interesting is the result under isotropic tension (fig. 8). Like in 
the previous cases we obtain a ramified main crack which is covered with 
many small side branches. Again we observe that theses side branches grow 
perpendicular to the main branch. This fact is particularly interesting in 
this geometry since here no growth direction is favored by the boundary 
condition and thus the right angles can only be explained by the previously 
given physical arguments. Furthermore one has to notice the stability of the 
crack tips which do not show a tendency to split. One rather observes the 
appearance of side branches behind the crack tip. 

4.3 Variation of the material strength 

In this section we want to discuss the effect of the material strength on the 
results of the simulation. We set the strength to two extreme values: either 
to zero or to ±100 while keeping M = oo fixed. 

In the figures 9-11 the results for a c — a t — are shown. Setting both 
values to zero means, that the crack grows as soon as it is touched by a 
walker. 

In the uniaxial tension experiment (fig. 9) the crack remains its elongated 
shape but the number and shape of the side branches changed drastically. 
They completely lost their straight appearance and are much more ramified. 
Also the pronounced 90° angles vanish. Now, the crack looks very much like 
a DLA cluster with anisotropic sticking probability . 

In the shear simulation (fig. 10) the global cross shape of the crack 
is still present. But here in contrast to the uniaxial case the right angles 
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between main and side branches and the very straight shape of all branches 
are still present. Amazingly, the dendritic shape of the cracks under shear is 
extremely stable and is not destroyed even in the limit of zero strength. 

In the case of isotropic tension (fig. 11) almost all features emphasized 
above are gone. The cluster does no longer reveal side branching and 90° 
angles between branches or screening of curved cracks. The resulting crack 
looks like a typical DLA cluster of the corresponding size. 

As a result one obtains, that in the limit of zero strength our model 
crosses over into a anisotropic DLA like behavior. This is readily explained 
because in the limit of zero strength the noise generated by the random 
walkers dominates over the structure imposed by the outer boundaries. But 
the vectorial character of the problem is not completely gone. It introduces 
an intrinsic anisotropy as is clearly seen by comparing the uniaxial and the 
shear case: In the first case one still obtains one straight crack perpendicular 
to the force, whereas in the latter case one obtains the cross like shear crack. 
Such a peculiarity cannot be found in the scalar DLA. 

In the figures 12-14 the results for a c — a t — 100 are shown. In the 
shear simulation (fig. 12) one observes that the dendritic shape of the crack 
is extremely stable. Even in the limit of high noise reduction the dendritic 
shape of the four main branches does not change. Neither the number nor 
the size of the main branches seems to decrease. Also in the case of uniaxial 
tension (fig. 13) the overall shape of the cluster does not change. One still 
observes a straight structure with many side branches and even the side 
branches that grow parallel to the main branch remain. In the isotropic 
geometry (fig. 14) one observes that the main and the side branches tend to 
be straighter and the overall crack seems to be less ramified. Here one can 
see clearly the tendency of side branches to grow perpendicular to the main 
branches. 

4.4 Asymmetric strength 

Here we want to address the question of asymmetric strengths. In typical 
real material like concrete one usually finds the behavior |o" c |>|cr t | . Thus, as 
special cases we want to comment on the situation, that the material only 
breaks under tension: a c = — oo, a t = 20 . We also studied the opposite case 
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in which the medium breaks only under compression: cr t = oo, a c = —20 . 

The most obvious change in the result is obtained in the shear experiments 
(fig. 15). As soon as one strength is chosen much larger than the other one, 
only one of the diagonal main branches remains. The curvature of the crack 
is purely a boundary effect. In fact, according to the breaking mode either 
one or the other diagonal is selected. This is in complete agreement with 
previous simulations . One also finds that the pronounced dendritic shape 
that was observed earlier has completely vanished and the single crack does 
no longer contain any side branch at all. 

A similar observation can be made in the uniaxial tension simulation (fig. 
16). If one only allows breaking under tension one obtains a single straight 
crack perpendicular to the force direction. Again, all side branches disappear. 
In the other case, in which only the breaking under compression is allowed, 
the crack perpendicular to the load direction disappears and a single crack 
parallel to the load direction is observed. This effect is readily understood 
with the Young experiment: A tensile load on an elongated bar will not only 
increase the length of the bar but it will also influence the width of the bar, 
according to the Poisson ratio. Thus, the breaking mode observed here is a 
result of the positive Poisson ratio. 

The suppression of compressive breaking seems to have the smallest effect 
in the isotropic case (fig. 17). Here, there seems to be no difference between 
a crack generated only for tensile breaking and a crack generated under 
tensile and compressive cracking: One observes 90° angles between main and 
side branches, one observes side branching rather than tip splitting and one 
obtains the screening behavior described earlier. 



4.5 Experiments 

We also reproduce with our new model real experiments performed on con- 
crete. As an example we want to show the results obtained for the "four- 
point-shear experiment" |32] and a mixed "tensile-shear experiment" |3l| . 



In the four point shear experiment a concrete sample containing two 
notches in it to initiate the fracture process at well defined points is loaded 
with four shear forces as indicated in figure 18. One has to mention, that in 
the experiments of course finite samples are used, which imposes additional 
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boundaries on the simulation. In our simulation we simplify this and apply 
four forces to an infinite plane. But we showed with a boundary-integral cal- 
culation of the stress field that this simplification does not lead to substantial 
difficulties. For simplification we describe the notches by two single parti- 
cles at the respective positions (fig. 18). The geometrical relations for the 
location and the magnitude of the forces are the same as in the experiment. 

The result of our simulation is shown in figure 19, whereas the result of 
experiments is shown in figure 20. Both results show the same behavior: One 
observes curved cracks, which start with an almost 45° angle and cross over 
into a straight behavior: In the initial stages of growth the crack develops 
according to a pure shear field generated by the two forces between which 
the microcracks are located. This results in an roughly 45° angle. If the 
cracks get larger the situation changes and growth occurs according to a 
compressive field generated by the two forces located near the notch. 

In the experiment one obtains in the end one large and one small crack, 
while we observe in the simulations a symmetric situation. One observes 
experimentally that initially two cracks start to grow, but one of them stops 



after one has reached the yield stress [32|. This yield stress is not accessible 



in our simulations, which is the reason for the symmetric crack pattern. 

In the second calculation we simulated another experiment also performed 
on a concrete sample which is known as load-path 2: To a concrete sample 
first a tensile load is applied until a crack of a certain size developed, then it 
is unloaded and a tensile shear load is applied. The result of our simulation 
(fig. 21) of this experiment is compared to the experimental results (fig. 22). 
One observes first a long straight tensile crack as a result of the tensile load. 
During the second load phase the diagonal shear cracks emerge. 



4.6 Biased walks 

In all simulations described above the walkers performed an isotropic random 
walk. This is in agreement with the calculations of Roux. In this section we 
want to consider the case of biased random walks, in which each walker jumps 
into a preferred direction. 

One physical motivation of the bias is the following. The stress g_ and the 
strain field e lead inside an elastic medium to an elastic energy U which is 
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stored in the fields ( 

U = ifiyVij ( 19 ) 

In our case of an infinite and isotropic elastic medium with a point force 

— * 

F at the origin this can be calculated at point f as 

U( r) oc 1 + 2 (2 ~* /)(1 t Z/) - cos 2 y? (20) 
(1 -v) 

where cos (p oc F-f. Now, we consider the regions with high elastic energy 
as hot regions in which crack growth is most likely. Thus it is reasonable to 
let the walker preferentially diffuse into these regions by using an angular 
jump probability of the type 

p((p) oc I + a- cos 2 (p (21) 

where a is now a free model parameter and ip is the angle between jump 
and force direction. 

Since we must not overestimate the contribution of each walker to the 
stress counters at the crack surface, we now have to accumulate scaled 
stresses 

where cosy? oc ( — 77) • F . Here we take into account the anisotropic 
distribution of walkers. 

In this formulation the walkers diffuse preferentially into the direction of 
the force, which is certainly a good choice for walkers that are launched at 
the outer boundary. But if one considers the walkers that are relaunched 
at the crack surface, this might lead to unphysical effects: Since all walkers 
diffuse until they touch the crack, the relaunched walkers might touch the 
crack again in an angle ip ~ n/2 which finally leads to a— >oo . Thus, we still 
overestimate the contribution of individual walkers. This shows, that the bias 
does not only increase the density of walkers in certain regimes of space, but 
it also shifts the balance between the walkers launched at the outer boundary 
and at the crack. Therefore, we also consider a second possibility of a bias, in 
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which the walkers preferentially diffuse perpendicular to the force direction 

p(ip) oc 1 + a- sin 2 ip (23) 



The figures 23 and 24 show typical crack structures generated with a 
cosine-type and sine-type bias under isotropic tension. Each structure was 
generated for a bias parameter a = 10, a Poisson ratio v = 0.2, material 
strength a ct = ±20 and a fatigue parameter M = 10. The structure for the 
cosine-bias contains 2, 000 particles whereas the structure for the sine-bias 
contains 6, 000 particles. 

The cosine-type bias leads to very regular needle like cracks (fig. 23) 
while the sine-type bias produces a completely new structure (fig. 24). 

Here, we observe very typical properties. The generated crack is very 
similar to the ones observed in the experiments by Lemaire and Van Damme 
With a sine-type bias one obtains a very pronounced 90° behavior: 



Initially all side branches grow perpendicular to their main branches. Like 
the cracks produced by Van Damme et al. the branches bend and side 
branches appear only on the outer side of the curved crack. Furthermore 
we observe, that again the crack tips are very stable and do not show a 
tip-splitting instability. 



4.7 Quantitative analysis 

All crack patterns shown above have a fractal structure which shall be ana- 
lyzed in terms of the dependence of the radius of gyration R g on the crack size 
N. Therefore we generated five clusters for each geometry uniaxial tension, 
shear and isotropic tension and each parameter set case 1.: high strength 
o ct = ±100 and strong fatigue M = 10,000 and case 2.: low strength 
o ct = ±30 and weak fatigue M = 10. The results are shown in the fig- 
ures 25-27. In all cases we find a fractal behavior 

N oc Rff (24) 

where the fractal dimensions Df are summarized in the following table. 

The fatigue-parameter M successfully suppresses side branching so that 
the resulting structures are essentially one dimensional as can be seen from 
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case 1. 


case 2. 


uniaxial tension 


1.89±0.02 


1.00±0.01 


shear 


1.51±0.01 


0.99±0.02 


isotropic tension 


1.56±0.02 


1.10±0.01 



the fractal dimensions in the second case: here the fractal dimensions are 
close to unity. In the case of strong fatigue — case 1. — the fractal dimen- 
sions for shear and isotropic tension are considerably lower than the fractal 
dimension of DLA, which is the result of the observed stability of the crack 
tips. Only the fractal dimension of the uniaxial-tension simulations is very 
large and indicates that the continuous side branching might lead to a com- 
pact structure. 

We also measured the fractal dimensions of the simulations using biased 
walks. For each case sine- and cosine-type bias we generated five cracks 
under isotropic tensile load. The results are shown in figure 28. Again, one 
observes a fractal behavior with dimensions 

„ J 1.03±0.05 cosine-type bias , . 

f ~ 1 1.23±0.04 sine-type bias ^ ' 

One observes that a cosine-type bias leads to one-dimensional, needle-like 
structures, whereas the sine-type bias leads to ramified fractals. 



5 Discussion 



As has been shown above, our new random walker method for crack growth 
successfully reproduces well known experimental results and results of previ- 
ous small scale simulations using lattice methods. 

The major advantage of our method is the complete absence of any lattice 
whatsoever. All three ingredients, the walk, the breaking criterion and the 
growth direction, have been formulated completely in a continuum language. 
Therefore all lattice effects — like the well known anisotropy in lattice-DLA 
35|j — are avoided. The observed anisotropics — like the pronounced cross- 



like structure under shear load — have a physical origin and are no artifacts. 
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This fact is particularly interesting for the observed right angles between 
main and side branches. We have to stress again that this off-lattice method 
reproduces this important physical property of real cracks. In fact, we are 
not aware of other simulations which reproduce this fact without using a 
structured lattice. 

An important feature of crack growth is the structural disorder of real 
materials. This disorder is implicitly included in our method by the use of 
random walkers, which introduces annealed disorder into the simulation. 

Another advantage of our method is that by using optimized algorithms 
[0] all calculations are easily done on workstations. Actually, our calculations 
are all performed on simple SPARC2 stations. They use about 5MB main 
memory for a 20,000 particle cluster. The typical run time for one cluster is in 
the order of 10-15 hours. By altering the strength and the fatigue parameter 
this time can be much shorter but also much longer. The interested reader 
can obtain all source codes from the author. 

One disadvantage of our method is that it is not easily possible to measure 
the displacement field u since the whole algorithm is based entirely on the 
stress field g_. Therefore it is not yet possible to measure stress-strain relations 
or energy release curves, which are of course important physical quantities. 
For the same reason it is also not possible to simulate experimental situations 
in which not the load but rather the displacement is controlled. However, 
one might argue that in the present stage one is restricted to infinitely stiff 
materials with vanishing displacements. 

The last point to be discussed is the fact that — except for the case of the 
biased walks — a walker undergoes an ordinary, unbiased random walk before 
touching the cluster. Thus, the probability that a walker reaches a certain site 
is implicitly governed by the Laplace equation. However, this does not impose 
a wrong, Laplacian screening behavior, which has two reasons. The first 
reason for this is, that the condition for relaunching a walker is completely 
determined by the Lame equation and the boundary condition and, thus, 
by the elasticity problem. Furthermore one has to notice that both, the 
Laplace and the Lame equation share the property, that the important fields, 
the electric and the stress field, diverge at tips. This is just an explanation of 
the fact that DLA clusters — like real cracks — preferentially grow at tips. 
This property is the reason why we chose random walkers to sample the hot 
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sites of the cracks. In our simulation we can even enhance this tip growth 
by strongly reducing the effect of fatigue and thereby inhibiting the growth 
of side branches. Thus, the fact that random walks are used to determine 
the "hot sites" of a crack is not in conflict with our aim to simulate crack 
growth. 

As an additional possibility we studied the case of biased walks, in which 
the bias is calculated from the applied forces. For this case we have generated 
cracks which are very similar to experimental ones observed by Lemaire and 



Van Damme |34]]. Although we do not yet fully understand the physical 
meaning of the bias parameter a we believe that further research into this 
direction will lead to new and interesting results. 



6 Summary 

In the present work we have studied the possibilities to implement a novel way 
to simulate crack growth by using random walkers. Inspired by Roux we have 
discussed stochastic algorithms — which are suggested by his original work — 
to solve the Lame equation. Although he showed an exact relation between 
random walkers and elasticity, a direct implementation of the method shown 
in his work is not practical, because it is not possible to suppress the strong 
fluctuations and to reach useful averages. However, his work showed that 
it is indeed possible to treat problems different from Laplacian ones with 
random walkers. Inspired by the closely related DLA algorithm — which 
has been very successful in solving problems of Laplacian growth — we have 
then constructed a different and new method to attack the crack growth in 
elastic media. 

Here, the random walkers are no longer used to construct a Green's func- 
tion, but are rather used to detect the hot sites in a crack, which are the next 
ones to grow. The formulation of the walk and the growth rule is purely made 
in the continuum to avoid all lattice effects. We are able to reproduce with 
this new method older results which have been obtained using other, conven- 
tional methods on smaller systems. We are also able to simulate experiments 
that have been performed on concrete. 

Although we restricted ourselves to the special case of Lame's equation in 
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the framework of linear elasticity, in principle our method can be extended 
to other growth phenomena. It may be interesting to investigate this aspect. 

I would like to thank H. J. Herrmann and S. Roux for many helpful 
discussions and comments. 



A Appendix 

Since the matrix M in the exponent of (H) is symmetric, it can be diagonal- 
ized and thus Q t { k) can be evaluated in closed form. One obtains in two 
dimensions 

\t( i\ _ 1 ( k 2 exp 2 +k 2 exp l A;i&;2{exp 2 — exp x } 



Q*( k) = "I;™ 1 " 2 ^\ 1,2 ,2 (26) 

= k 2 V kik 2 {exp 2 — exp x } k\ exp x +k 2 exp 2 ' 



(27) 



where exp^ = exp (— t\i k 2 ) and Ai and A2 are the eigenvalues 

2a + d 

1 = 2d(d + 2) 
3d+2a(l-d) 

2 ~ 2d(d + 2) 
of M . Now, the Fourier transformation into real space 

Q t (r) = [Q\k)exp(ikr)-^ (28) 
— J — (2tt) 

involves the calculation of integrals of the type 

J ^exp (-t\ k 2 ) exp (z k r)-0». (29) 

Using polar coordinates the evaluation of the angular integration leads to 
Bessel functions of the first and second kind, Jo(x) and J 2 (x) 

/■27T 

/ cos 2 if- exp (ik-r cos Up — ip r )) dip = ir[J (kr) — cos (2ip r ) J 2 (kr)] 
Jo 

r 

sin 2 if- exp (ik-r cos Up — <p r )) dp = it[Jo(kr) + cos (2ip r ) J 2 (kr)] (30) 

/ sin p> cos <f- exp (ik-r cos Up — <p r )) dip = tc[— sin (2ip r ) J 2 (kr)] 
Jo 
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The remaining radial integration involves an integral of Bessel functions 



and Gaussians and finally one obtains [36 



Q\ r) = hf s (r,t) + / B ( r ,t)].i - JL/ tt ( r ,0-2-^ (31) 
= Zn — Z7T r z 



where 



(32) 



exp (-r 2 /4Ait) exp (-r 2 /4A 2 t) 

/ 4AiA exp(-r 2 /4A 1 t) / 4A 2 A exp (-r 2 /4A 2 t) 
f a (r,t) = + ^ + 7^ 
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